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c3 Abstract 
Hi ■ 

An algorithm to compute polynomial conserved densities of polynomial nonlinear 
lattices is presented. The algorithm is implemented in Mathematica and can be used 
as an automated integrability test. With the code diffdens.m, conserved densities 
are obtained for several well-known lattice equations. For systems with parameters, 
the code allows one to determine the conditions on these parameters so that a 
' sequence of conservation laws exist. 
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^ ' 1 Introduction 

> 

There are several motives to find the explicit form of conserved densities of 
differential-difference equations (DDEs). The first few conservation laws have a 
physical meaning, such as conservation of momentum and energy. Additional 
ones facilitate the study of both quantitative and qualitative properties of 
solutions [1]. Furthermore, the existence of a sequence of conserved densities 
predicts integrability. Yet, the nonexistence of polynomial conserved quantities 
does not preclude integrability. Indeed, integrable DDEs could be disguised 
with a coordinate transformation in DDEs that no longer admit conserved 
densities of polynomial type [2]. 

Another compelling argument relates to the numerical solution of partial dif- 
ferential equations (PDEs). In solving integrable discretizations of PDEs, one 
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should check that conserved quantities indeed remain constant. In particular, 
the conservation of a positive definite quadratic quantity may prevent nonlin- 
ear instabilities in the numerical scheme. The use of conservation laws in PDE 
solvers is well-documented in the literature [3-5]. 



Several methods to test the integrability of DDEs and for solving them are 
available. Solution methods include symmetry reduction [6] and solving the 
spectral problem [7] on the lattice. Adaptations of the singularity confinement 
approach [8], the Wahlquist-Estabrook method [9], and symmetry techniques 
[10-12] allow one to investigate integrability. The most comprehensive integra- 
bility study of nonlinear DDEs was done by Yamilov and co-workers (see e.g. 
[2,13]). Their papers provide a classification of semi-discrete equations pos- 
sessing infinitely many local conservation laws. Using the formal symmetry 
approach, they derive the necessary and sufficient conditions for the existence 
of local conservation laws, and provide an algorithm to construct them. 

In [14,15], we gave an algorithm to compute conserved densities for systems 
of nonlinear evolution equations. The algorithm is based on the concept of 
invariance of equations under dilation symmetries. Therefore, it is inherently 
limited to polynomial densities of polynomial systems. Recently an extension 
of the algorithm towards DDEs was outlined in [16]. Here we present a full de- 
scription of our algorithm, which can be implemented in any computer algebra 
language. We also provide details about our software package diffdens.m [17] 
in Mathematica [18], which automates the tedious computation of closed-form 
expressions for conserved densities. 



Following basic definitions in Section 2, the algorithm is given in Section 3. 
To illustrate the method we use the Toda lattice [19], a parameterized Toda 
lattice [8], and the discretized nonlinear Schrodinger (NLS) equation [20-22]. 
In Section 4, we list results for an extended Lotka-Volterra equation [23], a 
discretized modified KdV equation [24], a network equation [24], and some 
other lattices [2]. The features, scope and limitations of our code diffdens.m 
are described in Section 5, together with instructions for the user. In Section 
6, we draw some conclusions. 



2 Definitions 
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2.1 Conservation laws 



Consider a system of DDEs which is continuous in time, and discretized in 
the single space variable, 

ii n = F(..., u n _i, u n , u n+ i, ...), (1) 

where u n and F are vector dynamical variables with any number of com- 
ponents. For simplicity of notation, the components of u„ will be denoted 
by u n ,v n , etc. We assume that F is polynomial with constant coefficients. If 
DDEs are non-polynomial or of higher order in t, we assume that they can be 
recast in the form (1). 

A local conservation law is defined by 

Pn — Jn ~ Jn+l, (2) 

which is satisfied on all solutions of (1). The function p n is the conserved 
density and J n is the associated flux. Both are assumed to be polynomials in 
u n and its shifts. 

Obviously, ^(E„Pn) = E„P« = Hn{Jn ~ J n +i), and this telescopic series 
vanishes for a bounded periodic lattice or a bounded lattice resting at infinity. 
In that case, J2n Pn is constant in time. So, we have a conserved quantity. 

Example 1 Consider the one-dimensional Toda lattice [19] 

y n = exp (y n _i - y n ) - exp (y n - y n+1 ), (3) 

where y n is the displacement from equilibrium of the nth unit mass under an 
exponential decaying interaction force between nearest neighbors. 

With the change of variables, 

u n = y n , v n = exp(y n -y n+1 ), (4) 

lattice (3) can be written in polynomial form 

Un = Vn-1 ~ V n , V n = V n (u n - U n+1 ). (5) 

System (5) is completely integrable [19,25]. The first two density-flux pairs 

pi 1} = u n , jjp = v n -i, and = \u\ + v n , J n 2) = u n v n _ 1: (6) 
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can be easily be computed by hand. 
2.2 Equivalence criterion 

We introduce a few concepts that will be used in our algorithm. Let D denote 
the shift-down operator and U the shift-up operator. Both are defined on 
the set of all monomials in u n and its shifts. If m is such a monomial then 
Dm = m| n ^ n _i and Um = m\ n ^ n+ i. For example, Du n+ 2V n = u n+ iv n -i and 
Uu n _ 2 v n -i — u n -iv n . It is easy to verify that compositions of D and U define 
an equivalence relation on monomials. Simply stated, all shifted monomials are 
equivalent, e.g. u n _iv n+ i = u n+2 v n+4: = -u„_ 3 t>„_i. This equivalence relation 
holds for any function of the dependent variables, but for the construction of 
conserved densities we will apply it only to monomials. 

In the algorithm we will use the following equivalence criterion: if two mono- 
mials mi and m 2 are equivalent, rxi\ = m 2 , then m\ = m 2 + [M n — M n+1 ] for 
some polynomial M n that depends on u„ and its shifts. For example, u n - 2 u n = 

Wn-l^n+l Since U n - 2 U n = U n ^iU n+1 + [u n ^ 2 U n - U n -iU n+ i] = M n _iM n+ i + [M n - 

M n+ i], with M n = u n _ 2 u n . 

The main representative of an equivalence class is the monomial of that class 
with n as lowest label on u (or v). For example, u n u n+2 is the main repre- 
sentative of the class with elements u n -iu n+ i, u n+ iu n+ 3, etc. Lexicographical 
ordering is used to resolve conflicts. For example, u n v n+2 (not u n _ 2 v n ) is the 
main representative in the class with elements u n _ 3 v n _i,u n+2 v n+i , etc. 



3 Algorithm 

Scaling invariance, which results from a special Lie-point symmetry, is an 
intrinsic property of many integrable nonlinear PDEs and DDEs. Indeed, ob- 
serve that (5), and the couples p£\ and in (6) (after inserting 
them in (2)), are invariant under the dilation (scaling) symmetry 

(t, u n , v n ) -> (At, X~ l u n , A~V), (7) 

where A is an arbitrary parameter. Stated differently, u n corresponds to one 
derivative with respect to t; denoted by u n ~ ^. Similarly, v n ~ 

Our three-step algorithm exploits the symmetry (7) to find conserved densities. 
Step 1: Determine the weights of variables 
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The weight, w, of a variable is by definition equal to the number of derivatives 
with respect to t the variable carries. Weights are non-negative, rational, and 
independent of n. Without loss of generality we set w(^) = 1. 

The rank of a monomial is defined as the total weight of the monomial. For 
now we assume that all the terms (monomials) in a particular equation have 
the same rank. This property is called uniformity in rank. The uniformity 
in rank condition leads to a system of equations for the unknown weights. 
Different equations in the vector equation (1) may have different ranks. 
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Example 2 Requiring uniformity in rank for each equation in (5) allows one 
to compute the weights of the dependent variables. Indeed, 

w(u n ) + 1 = w(v n ), w(v n ) + 1 = w(u n ) + w(v n ), (8) 



yields w(u n ) = 1, w(v n ) = 2, which is consistent with (7). 



Step 2: Construct the form of the density 

This step involves finding the building blocks (monomials) of a polynomial 
density with prescribed rank R. All terms in the density must have the same 
rank R. Since we may introduce parameters with weights (see Example 6), 
the fact that the density will be a sum of monomials of uniform rank does not 
necessarily imply that the density must be uniform in rank with respect to 
the dependent variables. 

Let V be the list of all the variables with positive weights, including parameters 
with weight. The following procedure is used to determine the form of the 
density of rank R: 

• Form the set Q of all monomials of rank R or less by taking all appropriate 
combinations of different powers of the variables in V. 

• For each monomial in Q, introduce the appropriate number of derivatives 
with respect to t so that all the monomials exactly have weight R. Gather 
in set Ti all the terms that result from computing the various derivatives. 

• Identify the monomials that belong to the same equivalence classes and 
replace them by the main representatives. Call the resulting simplified set 
X, which consists of the building blocks of the density with desired rank R. 

• Linear combination of the elements in X with constant coefficients q gives 
the form of polynomial density of rank R. 

Example 3 Continuing with (5), we compute the form of the density of rank 
R = 3. From V = {u n , v n } we build Q = {u n , u n , u n v n , u n , v n }. Next, introduce 
^-derivatives, so that each term exactly has rank 3. Thus, using (5), 



d° d° 

/ 3 \ 3 i \ 

~^0\ U n) = \) ~^oV U n v n) = U n V n , 

^( U n 2 ) = 2u n U n = 2u n V n - 1 - 2u n V n , "j"^™) = *»» = ~ U n+1 V n , 



d 

dt 2y ~ ,l/ dv~ ,v dt 



2 \ U n) = — («n) = — (Vn-1 - V n ) = Mn-l^n-1 ~ U n V n ^ - U n V n + U n+1 V n . 



Gather the resulting terms in set Ti = {u n , u n v n -i, u n v n , w„_iv n _i, u n+ iv n }. 
Since u n +iv n = u n v n -\ and w n _if n _i = u n v n in H, u n+ iv n is replaced by 
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u n v n -i, and u n ~iv n -i by u n v n . Hence, we obtain J = {u n , u n v n -i, u n v n }. 
Linear combination of the monomials in I gives the form of the density: 

Pn = Ci u\ + C 2 M„V„-l + C 3 U n V n . (9) 

S'tep 5/ Determine the unknown coefficients in the density 

The following procedure simultaneously determines the constants q and the 
form of the flux J n : 

• Compute p n and use (1) to remove all the t-derivatives. 

• Regarding (2), the resulting expression must match the pattern J n — J n +i- 
Use the equivalence criterion to modify p n . The goal is to introduce the 
main representatives and to identify the terms that match the pattern. 

• The part that does not match the pattern must vanish identically for any 
combination of the components of u n and their shifts. This leads to a linear 
system S in the unknowns q. If S has parameters, careful analysis leads to 
conditions on these parameters guaranteeing the existence of densities. See 
[14] for a description of this compatibility analysis. 

• The flux is the first piece in the pattern [J n — J n+ i] . 

Example 4 Carrying on with (5), we determine the coefficients C\ through 
c 3 in (9) by requiring that (2) holds. Simultaneously, we determine J n . 

Compute p n using (9). Use (5) to remove u n , v n , etc. After grouping the terms 

pn = (3Ci - C 2 )u 2 n V n -l + (C 3 - Sc^ulVn + (c 3 - C 2 )V n _ 1 V n 

+c 2 u n ^ 1 u n v n ^i + c 2 v 2 n _ 1 - c 3 u n u n+1 v n - c 3 vl. (10) 

Use the equivalence criterion to modify p n . For instance, replace u n -\U n v n -\ 
by u n u n+ iv n + [u ri -iu n v n -i — u n u n+ iv n ]. In terms of main representatives, 

p n = (3ci - c 2 )u n v n -\ + (ca - 3ci)«^u n 

+ (C 3 - C 2 )v n V n+ i + [(C 3 - C 2 )v n „iV n - (C 3 ~ C 2 )v n Vn+l] 
+C 2 U n U n+1 V n + [c 2 u 

n—~lV"nVn—l C 2 U n Un+\V n \ 
+C 2 V 2 n + [C 2 ^_! - C 2 vl\ - C 3 U n U n+1 V n - C 3 V^. (11) 

Next, group the terms outside of the square brackets and move the pairs inside 
the square brackets to the bottom. Rearrange the latter terms so that they 
match the pattern [J n — J n +i\- Hence, 

pn = (3Ci - C 2 )u\v n ~\ + (c 3 - ?>C X )u n V n 
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+ (c 3 - c 2 )v n v n+1 + (c 2 - c 3 )u n u n+1 v n + (c 2 - c 3 )^ 

+ [{(C 3 - C 2 )v n _iV n + C 2 M n _iW„V„_i + C 2 vl_ 1 } 

-{(c 3 - c 2 )w n Vn+i + C 2 M n M„ + it;„ + C 2 ^}]. (12) 

The first piece inside the square brackets determines 

</n = (C 3 - C 2 )v n - X Vn + C 2 U n - 1 U n V n „ 1 + C 2 vl_ v (13) 

The terms outside the square brackets must all vanish piece by piece, yielding 
S = {3 Cl - c 2 = 0, c 3 - 3ci = 0, c 2 - c 3 = 0}. (14) 

The solution is 3ci = c 2 = c 3 . Since densities can only be determined up to 
a multiplicative constant, we choose C\ = |, c 2 = c 3 = 1, and substitute this 
into (9) and (13). The explicit forms of the density and the flux follow: 

Pn = I ul + U n (v n ^ + V n ), J n = U n - X U n V n - X + V 2 n _ v (15) 

Example 5 To illustrate how the algorithm works for DDEs with parameters, 
consider the parameterized Toda lattice 

Un = a Vn-! - v n , v n = v n ((3 u n - u n+1 ) , (16) 

where a and (3 are nonzero parameters without weight. In [8] it was shown 
that (16) is completely integrable if and only if a — f3 — 1; but then (16) is 
(5). 

Using our algorithm, one can easily compute the compatibility conditions for 
a and /3, so that (16) admits a polynomial conserved densities of, say, rank 3. 
The steps are the same as for (5). However, (14) must be replaced by 

S = {3acx — c 2 = 0, /3c 3 — 3c x = 0, ac 3 — c 2 = 0, (3c 2 — c 3 = 0, ac 2 — c 3 = 0}. 



A non-trivial solution 3ci = c 2 = c 3 will exist if and only if a — f3 — 1. 

Analogously, (16) has density p n l > = u n of rank 1 if a — 1, and density 
P^n = 2 u n + v n of rank 2 if af3 = 1. Only when a — — 1 will (16) have 
conserved densities of rank > 3: 



p^ 3) = \u n + w„(v„_i + v n ), 

P ( n = Wn + + V n ) + U n U n+l V n + W n + U n V n+ i, 



(17) 

(18) 



8 



/4 5) = \u b n + v? n {v n -i + v n ) + u n u n+1 v n (u n + U n+ i) 

+U n V n _ 1 (v n _ 2 + + ^n) + WA^n-l + + «n+l)- (19) 

Ignoring irrelevant shifts in n, these densities agree with the results in [25]. 

Example 6 In [20,21], Ablowitz and Ladik studied the following integrable 
discretization of the NLS equation: 

i u n = u n+1 - 2u n + M n _i + <« n (« n+ i + lin_i), (20) 

where it* is the complex conjugate of u n . Instead of splitting u n into its real 
and imaginary parts, we treat u n and v n = u* n as independent variables and 
augment (20) with its complex conjugate. Absorbing % in the scale on t, 

u n = u n+l ~ 

2u n + U n _i + u n v n (u n+1 + M„_i), 
Vn = -(V n +1 ~ 2V„ + V n -i) ~ U n V n {v n +l + U n -l)- (21) 

Since v n = u* n we have w(v n ) = w(u n ). Neither of the equations in (21) 
is uniform in rank. To circumvent this problem we introduce an auxiliary 
parameter a with weight, and replace (21) by 

Un = «(«n+l - 2«„ + + U n V n (u n+1 + 

u n = -a(v n+1 - 2v n + w„_i) - u n v n (v n+1 + v n _i). (22) 
This extra freedom allows us to impose uniformity in rank: 

w(u n ) + 1 = w(a) + w(u n ) = 2w(u n ) + w(v n ) = 3w(u n ), (23) 
w(v n ) + 1 = w(a) + w(v n ) = 2w(v n ) + w(u n ) = 3w(v n ), (24) 

which yields w(u n ) = w(v n ) = |,u>(ct) = 1, or, u 2 n ~ v 2 ~ a ~ ^. 

We show how to get the building blocks of the density of rank |. In this case 
V = {u n ,v nj a} and = {u n ,v nj a,u 2 n ,u n v n ,vl,au n ,u 3 n ,av n ,u 2 n v n ,u n vl,v 3 }. 
The monomials cm n , w^, av n , u 2 n v n , u n v 2 and have already rank |, so no 
derivatives are needed. The monomials u n and v n will have rank | after in- 
troducing ^j. There is no way for the remaining monomials a, w n 2 , and 
f„ 2 to have rank | after differentiation with respect to t. Therefore, they are 
rejected. The remaining intermediate steps lead to 

J = {au n , ul, av n , u 2 n v n , U n U n+ iV n , u n v n - X v n , u n v 2 n , v*, u n u n+1 v n+1 , U n V n V n+ i}. 

Although uniformity in rank is essential for the first two steps of the algorithm, 
after the second step, we may already set a — 1. The computations now 
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proceed as in the previous example. The trick of introducing one or more 
extra parameters with weights can always be attempted if any equation in (1) 
lacks uniformity in rank. 

We list some conserved densities of (21): 

pil ] = C\U n v n -\ + c 2 u n v n+1 , (25) 

P ( n = c l{\ u l v l-l + U n U n+1 V n - X V n + U n V n ^ 2 ) 

+ c 2 {\u 2 n v 2 n+1 + u n u n+1 v n+1 v n+2 + u n v n+2 ), (26) 

Pn ] = C lil U l V l-l + MnM„+lV„-lV„(M re V„-l + U n+l V n + U n+2 V n+1 ) 

+ U n V n -i(u n V n ^ 2 + U n+ iV n -i) + U n V n (u n+ iV n - 2 + U n+2 V n ^i) + U n V n - 3 ] 
+ c 2[|«n^n+l + U n U n+l V n+l V n+2 {u n V n+ i + U n+1 V n+2 + U n+2 V n+3 ) 

+u n v n+2 (u n v n+l + u n+1 v n+2 )+u n v n+3 (u n+1 v n+1 +u n+2 v n+2 )+u n v n+3 \. (27) 

Our results confirm those in [20]. Also, if defined on an infinite interval, (20) 
admits infinitely many independent conserved densities [20]. Although it is a 
constant of motion, we cannot find the Hamiltonian of (20), 

H = -i EKK-i + u n+1 ) - 2 ln(l + u n u* n )l (28) 

n 

for it has a logarithmic term [22]. 



4 More examples 

4-1 An extended Lotka-Volterra equation 

In [23] , Hu and Bullough considered an extended version of the Lotka-Volterra 
equation: 

fe-i 

«n = J2( U n-r ~ U n+r )u n . (29) 
r=l 

For k = 2, (29) is the well-known Lotka-Volterra equation, for which the 
densities were presented in [16]. 

We computed five densities of (29) for the case k = 3. The first three are: 
Pi = u n , p2 = \u 2 n + u n (u n+l + u n+2 ), (30) 
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P3 = ^U 3 n + U 2 n {u n+1 + U n+2 ) + U n {u 2 n+1 + U 2 n+2 ) 

+U n U n+1 (u n+2 + U n+3 ) + U n U n+2 (u n+ i + U n+3 + U n+4 ). (31) 

For k = 4, we also computed five densities of (29). The first three are: 



Pi = Un, p 2 = ^U 2 n + U n (u n+1 + U n+2 + U n+3 ) , (32) 

P3 = \u\ + u l(u n+1 + U n+2 + U n+3 ) + U n (u n+ i + U n+2 + U n+3 ) 2 

+ u n u n+3 (u n+4 + u n+5 + u n+6 ) + u n u n+A {u n+1 + u n+2 ) + u n u n+2 u n+5 . (33) 



We computed four densities of (29) for k = 5. To save space we do not list 
them here. The integrability and other properties of (29) are discussed in [23] . 



4-2 A discretized modified KdV equation 

In [24], we found the following integrable discretization of the MKdV equation: 

tin = (1 + u n )(u n+1 - m„_i). (34) 

We computed four densities of (34). The first three are: 

= u n u n+1 , p {2) = \u 2 n u 2 n+l + u n u n+2 {l + u 2 n+1 ), (35) 

Pn ] = \ U l U l+l + U n U n+1 U n+2 (u n + U n+2 ){1 + U 2 n+l ) 

U n+2)- 

(36) 

4-3 Self-dual network equations 

The integrable nonlinear self-dual network equations [2,24] can be written as: 

ii n = (1 + u 2 n )(v n - u„_i), v n = {l + v 2 n ){u n+ i-u n ). (37) 

We computed the first four densities of (37). The first three are 

=u n v n _ 1 +u n v n , (38) 

P { n = \ U liyl-l + V l) + U n U n+l (l + V 2 n ) + V n {u 2 n V n _ X + V n+1 ) , (39) 
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= \ul{vl_ x + vl) + U n U n+ i{u n V n -i + M„ + iW„ + u n v n )(l + U*) 

+u n v n+ i(l + «^ +1 )(l + 0- ( 4 °) 
^.^ Generalized lattices 

Shabat and Yamilov [2] studied the following integrable Volterra lattice: 

il n = U n (v n+ i - V n ) , V n = V n (u n - Un_i). (41) 

With our program we computed the first four densities for this system: 

Pn )= «n + Un, p£ 2) = I («n + V n) + "n( v n + , (42) 

p^ 3) = |(w 3 l + t' 3 )+^(v n + v n+ i)+w n (^ + ^ +1 ) + w n i' n+ i(w„ + i + v n ), (43) 

Pi 4) = l(< + O + «nK + V n+1 ) + |<K + 

+u n (v n + v n+1 ) + 2u n v n+1 (u n + u n+l) 

+u n u n+ iv n+1 (u n + u n+ i +v n + v n+2 ) + u n v n v n+ i(v n + v n+ i). (44) 
In [2], the following Hamiltonian lattice is also listed: 

U n = U n+ i + U 2 n V n , V n = -V n -i - U n V n . (45) 

Four densities of (45) are: 

p n 1] = u n v n , p {2) = \u 2 n v 2 n + U n V n ^, (46) 

p^> = \u n v n + u n v n (u n v n -i + U n+1 V n ) + U n V n - 2 , (47) 

Pn ) = \ U n V n + ^(Wl + Mn+l^n) + + «n^n-3 

+u n v n (u n v n - 2 + 2u n+ iv n ^i + M n+2 v„ + i4 +1 u„i;„+i). (48) 
5 The Mathematica code diffdens.m 

We describe the features, scope and limitations of our program diffdens.m, 
which is written in Mathematica syntax [18]. The program has its own menu 
interface with a dozen data files. Users should have access to Mathematica. 
The code diffdens.m and the data files [17] must be in the same directory. 
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5.1 The menu interface 

After launching Mathematica, type 
In[l] := «diffdens.m 

to read in the code diffdens.m and start the program. Via its menu interface, 
the program will prompt the user for answers. 

The density is available at the end of the computations. To view it in standard 
Mathematica notation, type rho. To display it in a more elegant subscript 
notation, type subscriptf orm[rho] . 

5. 2 Preparing data files 

To test systems that are not in the menu, one should prepare a data file in 
the format of the data files that are provided with the code. Of course, the 
name for a new data file should not coincide with any name already listed in 
the menu, unless one intended to modify those data files. 
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Example 7 For the parameterized Toda lattice (16) the data file reads: 

(* Start of data file d_ptoda.m with parameters. *) 
(* Toda Lattice with parameters aa and bb *) 

u[l] [n_] ' [t] := aa*u[2] [n-1] [t]-u[2] [n] [t] ; 

u[2] [n_] ' [t] := u[2] [n] [t]*(bb*u[l] [n] [t]-u[l] [n+1] [t] ) ; 

noeqs = 2; 

name = "Toda Lattice (parameterized)"; 
parameters = -[aa,bb}; 
weightpars = O; 

(* The user may give the rank of rho *) 
(* and a name for the output file. *) 
(* rhorank = 3; *) 
(* myfile = "ptodar3.o"; *) 

(* The user can give the weights of u[l] and u[2] , *) 
(* and of the parameters with weight (if applicable) . *) 
(* weightu[l] = 1; weightu[2] = 2; weight [aa] = 1; *) 

formrho = 0; 

(* The user can give the form of rho. *) 
(* For example, for the density of rank 3: *) 
(* formrho = c[l]*u[l] [n] [t] ~3+c [2] *u [1] [n] [t]*u[2] [n-1] [t] + 
c[3]*u[l] [n] [t]*u[2] [n] [t] ; *) 

(* End of data file d_ptoda.m *) 

A brief explanation of the lines in the data file now follows. 
u[i][n_]'[t] := 

Give the i th equation of the system in Mathematica notation, 
noeqs = 2; 

Specify the number of equations in the system, 
name = "Toda Lattice (parameterized)"; 

Pick a short name for the system under investigation. The quotes are essential. 
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parameters = {aa,bb}; 

Give the list of parameters in the system. If none, set parameters = { }. 
weightpars = O; 

Give the list of the parameters that are assumed to have weights. Note that 
weighted parameters are not listed in parameters, the latter is the list of pa- 
rameters without weight. 

(* rhorank =3; *) 

Optional. Give the desired rank of the density, if less interactive use of the 
program is preferred (batch mode). 

O myfile = "ptodar3.o"; *) 

Optional. Give a name of the output file, again to bypass interaction with the 
program. 

O weightu[l] = 1; weightu[2] = 2; *) 

Optional. Specify a choice for some or all of the weights. The program then 
skips the computation of the weights, but still checks for consistency. Particu- 
larly useful if there are several free weights and non-interactive use is preferred. 

formrho = 0; 

If formrho is set to zero, the program will compute the form of p n . 

formrho = c[l]*u[l] [n] [t] ~3+c [2] *u [1] [n] [t]*u[2] [n-1] [t] + 
c[3]*u[l] [n] [t]*u[2] [n] [t] ; 

Alternatively, one could give a form of p n (here of rank 3). The density must 
be given in expanded form and with coefficients c[i]. If form of p n is given, 
the program skips both the computation of the weights and the form of the 
density. Instead, the code uses what is given and computes the coefficients c[i]. 
This option allows one, for example, to test densities from the literature. 

Anything within (* and *) is a comment and ignored by Mathematica. 

Once the data file is ready, pick the choice "tt) Your System" in the menu. 
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5.3 Scope and limitations 



Our program can handle systems of first order DDEs that are polynomial 
in the dependent variables. Only one independent variable (£) is allowed. No 
terms in the DDEs should have coefficients that explicitly depend on t or n. 
The program only computes polynomial conserved densities in the dependent 
variables and their shifts, without explicit dependencies on t or n. 

Theoretically, there is no limit on the number of DDEs. In practice, for large 
systems, the computations may take a long time or require a lot of memory. 
The computational speed depends primarily on the amount of memory. 

By design, the program constructs only densities that are uniform in rank. 
The uniform rank assumption for the monomials in p n allows one to compute 
independent conserved densities piece by piece, without having to split lin- 
ear combinations of conserved densities. Due to the superposition principle, a 
linear combination of conserved densities of unequal rank is still a conserved 
density. This situation arises frequently when parameters with weight are in- 
troduced in the DDEs. 

The input systems can have one or more parameters, which are assumed to be 
nonzero. If a system has parameters, the program will attempt to compute the 
compatibility conditions for these parameters such that densities (of a given 
rank) exist. The assumption that all parameters in the given DDE must be 
nonzero is essential. As a result of setting parameters to zero in a given system 
of DDEs, the weights and the rank of p n might change. 

In general, the compatibility conditions for the parameters could be highly 
nonlinear, and there is no general algorithm to solve them. The program au- 
tomatically generates the compatibility conditions, and solves them for pa- 
rameters that occur linearly. Grobner basis techniques could be used to re- 
duce complicated nonlinear systems into equivalent, yet simpler, non-linear 
systems. For DDEs with parameters and when the linear system for the un- 
known coefficients q has many equations, the program saves that system and 
its coefficient matrix, etc., in the file worklog.m. Independent from the pro- 
gram, the worklog files can later be analyzed with appropriate Mathematica 
functions. 

The assumption that the DDEs are uniform in rank is not very restrictive. 
If the uniform rank condition is violated, the user can introduce one or more 
parameters with weights. This also allows for some flexibility in the form of 
the densities. Although built up with terms that are uniform in rank in the 
dependent variables and parameters, the densities do no longer have to be uni- 
form in rank with respect to the dependent variables alone. Conversely, when 
the uniform rank condition is fulfilled, the introduction of extra parameters 
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(with weights) in the given DDE leads to a linear combination of conservation 
laws, not to new ones. 

In cases where it is not clear whether or not parameters with weight should be 
introduced, one should start with parameters without weight. If this causes 
incompatibilities in the assignment of weights (due to non-uniformity), the 
program may provide a suggestion. Quite often, it recommends that one or 
more parameters be moved from the list parameters into the list weightpars of 
weighted parameters. 

For systems with two or more free weights, the user will be prompted to 
enter values for the free weights. If only one weight is free, the program will 
automatically compute some choices for the free weight, and continue with the 
lowest integer or fractional value. The program selects this value for the free 
weight; it is just one choice out of possibly infinitely many. If the algorithm 
fails to determine a suitable value, the user will be prompted to enter a value 
for the free weight. 

Negative weights are not allowed. Zero weights are allowed, but at least one 
of the dependent variables must have positive weight. The code checks this 
requirement, and if it is violated the computations are aborted. Note that 
fractional weights and densities of fractional rank are permitted. 

Our program is a tool in the search of the first half-dozen conservation laws. 
An existence proof (showing that there are indeed infinitely many conservation 
laws) must be done analytically. If our program succeeds finding a large set 
of independent conservation laws, there is a good chance that the system of 
DDEs has infinitely many conserved densities. If the number of conserved 
densities is 3 or less, the DDE may have other than polynomial conserved 
densities, or may not be integrable (in the chosen coordinate representation). 



6 Conclusions 

We offer the scientific community a Mathematica package to carry out the 
tedious calculations of conserved densities for systems of nonlinear DDEs. 
The code diffdens.m, together with several data and output files, is available 
via Internet URL [17]. 

For lattices with parameters, the code automatically determines the compat- 
ibility conditions on these parameters so that a sequence of polynomial con- 
served densities exists. 

The existence of a large number of conservation laws is an indicator of inte- 
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grability of the lattice. Therefore, by generating the compatibility conditions, 
one can analyze classes of parameterized DDEs and filter out the candidates 
for complete integrability. 

Future generalizations of the algorithm will exploit other symmetries in the 
hope to find conserved densities of non-polynomial form. 

Acknowledgements 

Unal G6kta§ thanks Wolfram Research, Inc. for an internship. We acknowledge 
helpful comments from Drs. B. Herbst, S. Mikhailov, W.-H. Steeb, Y. Suris, 
P. Winter nitz, and R. Yamilov. We also thank G. Erdmann for his help with 
this project. 

References 

[I] E.G.B. Hohler and K. Olaussen, Int. J. Mod. Phys. A 11 (1996) 1831. 
[2] A.B. Shabat and R.I. Yamilov, Leningrad Math. J. 2 (1991) 377. 

[3] F.J. Hickernell, Stud. Appl. Math. 69 (1983) 23. 

[4] R.J. LeVeque, Numerical methods for conservation laws (Birkhauser Verlag, 
Basel, 1992). 

[5] J.M. Sanz-Serna, J. Comput. Phys. 47 (1982) 199. 

[6] D. Levi and P. Winternitz, J. Math. Phys. 34 (1993) 3713. 

[7] D. Levi and O. Ragnisco, Lett. Nuovo Cimento 22 (1978) 691. 

[8] A. Ramani, B. Grammaticos and K.M. Tamizhmani, J. Phys. A: Math. Gen. 
25 (1992) L883. 

[9] B. Deconinck, Phys. Lett. A 223 (1996) 45. 

[10] I. Cherdantsev and R. Yamilov, in: Symmetries and integrability of difference 
equations, eds. D. Levi, L. Vinet and P. Winternitz (American Mathematical 
Society, Providence, Rhode Island, 1996) 51. 

[II] D. Levi and R. Yamilov, J. Math. Phys. 38 (1997) 6648. 

[12] V.V. Sokolov and A.B. Shabat, in: Soviet Sci. Rev. Sec. C, Vol. 4 (Harwood 
Academic Publishers, New York, 1984) 221. 

[13] R. Yamilov, in: Proc. 8th int. workshop on nonlinear evolution equations and 
dynamical systems (World Scientific Publishing, Singapore, 1993) 423. 



18 



[14] U. Goktas, and W. Hereman, J. Symb. Comp. 24 (1997) 591. 

[15] U. G6kta§ and W. Hereman. The program condens.m is available via Internet 



URL: bttp : / / www . mines . edu/ fs_home / wher eman / 



[16] U. G6kta§, W. Hereman and G. Erdmann, Phys. Lett. A 236 (1997) 30. 

[17] U. G6kta§ and W. Hereman. The program diffdens.m is available from the 
WWW site in [15]. 

[18] S. Wolfram, The Mathematica book. 3rd Edition (Wolfram Media, Urbana- 
Champaign, Illinois & Cambridge University Press, London, 1996). 

[19] M. Toda, Theory of nonlinear lattices (Springer Verlag, Berlin, 1981). 

[20] M.J. Ablowitz and J.F. Ladik, J. Math. Phys. 17 (1976) 1011. 

[21] M.J. Ablowitz and J.F. Ladik, Stud. Appl. Math. 55 (1976) 213. 

[22] M.J. Ablowitz and B.M. Herbst, SIAM J. Appl. Math. 50 (1990) 339. 

[23] X-B. Hu and R.K. Bullough, J. Phys. A: Math. Gen. 30 (1997) 3635. 

[24] M.J. Ablowitz and P.A. Clarkson, Solitons, nonlinear evolution equations and 
inverse scattering (Cambridge University Press, Cambridge, U.K., 1991). 

[25] M. Henon, Phys. Rev. B 9 (1974) 1921. 



19 



